function J = f_jacobian(fun, theta)
    norm = sqrt(sum(theta.^2));
    tol  = max(norm, 1) * eps^(1 / 3);
    nk   = length(theta);
    J    = zeros(length(fun(theta)), nk);
    for i = 1:nk
        thetal     = theta; thetar = theta;
        thetal(i)  = thetal(i) - tol;
        thetar(i)  = thetar(i) + tol;
        J(:, i)    = (fun(thetar) - fun(thetal)) / (2 * tol);
    end
end
